function [V1, T1, E1, TE1, ET1] = refine_mesh(V, T, E, TE, ET, refine_E)
% function [V1, T1, E1, TE1, ET1] = refine_mesh(V, T, E, TE, ET, refine_E)
%
% To be finished with the return value, not globally!!
% This is the serial version.
%
%
%

% step 1: generate additional nodes:
% this is a temp flag for each old edge
idx_E = find(refine_E==1); nv = length(idx_E);
% add new nodes for each edge to be refined
idx_V = ((posV+1):(posV+nv))'; posV = posV + nv;
V(idx_V,:) = (V(E(idx_E,1),:)+V(E(idx_E,2),:))/2;
% fill the flag for E
flag_E = zeros(posE,3); flag_E(idx_E,3) = idx_V;

% step 2: seperate refined triangle into three categories:
% calculate the number of edge wait to be refined for each triangle
flag_T = refine_E(TE(1:posT,1)) + refine_E(TE(1:posT,2)) + refine_E(TE(1:posT,3));
% get the three categories of triangle:
T1 = find(flag_T==1); T2 = find(flag_T==2); T3 = find(flag_T==3);

% step 3: treate in three categories:
% run over the first category:
for t = 1:length(T1)
    tri = T1(t);
    eg_idx = find(refine_E(TE(tri,:))==1); % must be one element
    loc1 = eg_idx; loc2 = mod(loc1,3)+1; loc3 = mod(loc2,3)+1;
    v1 = T(tri,loc1); v2 = T(tri,loc2); v3 = T(tri,loc3);
    % A: add triangle first
    tri1 = posT+1; tri2 = tri; posT = posT+1; 
    % B: cope with edge
    eg = TE(tri,eg_idx); vc = flag_E(eg,3); % eg is global index
    if flag_E(eg,1) == 0 % this edge have not been divided yet, refine it!
        eg1 = posE+1; eg2 = eg; posE = posE + 1;
        E(eg1,:) = [v2,vc];  E(eg2,:) = [v3,vc];
        ET(eg1,1) = tri1; ET(eg2,1) = tri2;
        flag_E(eg,1) = eg1;  flag_E(eg,2) = eg2; 
    else % this edge has been divided already
        eg1 = flag_E(eg,2);  eg2 = flag_E(eg,1);
        ET(eg1,2) = tri1; ET(eg2,2) = tri2;
    end
    eg3 = posE + 1; posE = posE + 1; E(eg3,:) = [v1,vc];
    ET(eg3,:) = [tri1,tri2];
    e1 = TE(tri,loc1); e2 = TE(tri,loc2); e3 = TE(tri,loc3);
    ET(e2,find(ET(e2,:)==tri)) = tri2;
    ET(e3,find(ET(e3,:)==tri)) = tri1;
    % C: refine this triangle
    T(tri1,:) = [v1,v2,vc];      T(tri2,:) = [v1,vc,v3];
    TE(tri1,:) = [eg1,eg3,e3];TE(tri2,:) = [eg2,e2,eg3];
end

% run over the second category:
for t = 1:length(T2)
    tri = T2(t);
    eg_idx = find(refine_E(TE(tri,:))==0); % must be one element
    loc1 = eg_idx; loc2 = mod(loc1,3)+1; loc3 = mod(loc2,3)+1;
    v1 = T(tri,loc1); v2 = T(tri,loc2); v3 = T(tri,loc3);
    e1 = TE(tri,loc1); e2 = TE(tri,loc2); e3 = TE(tri,loc3);
    len2 = (V(v1,1)-V(v3,1))^2 + (V(v1,2)-V(v3,2))^2;
    len3 = (V(v1,1)-V(v2,1))^2 + (V(v1,2)-V(v2,2))^2;
    % A: add triangle first
    tri1 = posT+1; tri2 = posT+2; tri3 = tri; posT = posT+2; 
    if len2 >= len3
        % B: cope with edge
        eg = zeros(6,1);  
        if flag_E(e2,1) == 0 % this edge have not been divided yet, refine it!
            eg(3) = posE+1; eg(4) = e2; posE = posE + 1;
            E(eg(3),:) = [v3,flag_E(e2,3)];  E(eg(4),:) = [v1,flag_E(e2,3)];
            ET(eg(3),1) = tri3; ET(eg(4),1) = tri1;
            flag_E(e2,1) = eg(3);  flag_E(e2,2) = eg(4); 
        else % this edge has been divided already
            eg(3) = flag_E(e2,2);  eg(4) = flag_E(e2,1);
            ET(eg(3),2) = tri3; ET(eg(4),2) = tri1;
        end
        if flag_E(e3,1) == 0 % this edge have not been divided yet, refine it!
            eg(5) = posE+1; eg(6) = e3; posE = posE + 1;
            E(eg(5),:) = [v1,flag_E(e3,3)];  E(eg(6),:) = [v2,flag_E(e3,3)];
            ET(eg(5),1) = tri1; ET(eg(6),1) = tri2;
            flag_E(e3,1) = eg(5);  flag_E(e3,2) = eg(6); 
        else % this edge has been divided already
            eg(5) = flag_E(e3,2);  eg(6) = flag_E(e3,1);
            ET(eg(5),2) = tri1; ET(eg(6),2) = tri2;
        end
        eg(1) = posE + 1; eg(2) = posE + 2; posE = posE + 2;
        E(eg(1),:) = [flag_E(e2,3),flag_E(e3,3)];
        E(eg(2),:) = [flag_E(e2,3),v2];
        ET(eg(1),:) = [tri1;tri2]; ET(eg(2),:) = [tri2;tri3];
        ET(e1,find(ET(e1,:)==tri)) = tri3;
        % C: refine this triangle
        T(tri1,:) = [v1,flag_E(e3,3),flag_E(e2,3)];  TE(tri1,:) = [eg(1),eg(4),eg(5)];
        T(tri2,:) = [v2,flag_E(e2,3),flag_E(e3,3)];  TE(tri2,:) = [eg(1),eg(6),eg(2)];
        T(tri3,:) = [v3,flag_E(e2,3),v2];  TE(tri3,:) = [eg(2),e1,eg(3)];
    else
        % B: cope with edge
        eg = zeros(6,1);  
        if flag_E(e2,1) == 0 % this edge have not been divided yet, refine it!
            eg(3) = posE+1; eg(4) = e2; posE = posE + 1;
            E(eg(3),:) = [v3,flag_E(e2,3)];  E(eg(4),:) = [v1,flag_E(e2,3)];
            ET(eg(3),1) = tri2; ET(eg(4),1) = tri1;
            flag_E(e2,1) = eg(3);  flag_E(e2,2) = eg(4); 
        else % this edge has been divided already
            eg(3) = flag_E(e2,2);  eg(4) = flag_E(e2,1);
            ET(eg(3),2) = tri2; ET(eg(4),2) = tri1;
        end
        if flag_E(e3,1) == 0 % this edge have not been divided yet, refine it!
            eg(5) = posE+1; eg(6) = e3; posE = posE + 1;
            E(eg(5),:) = [v1,flag_E(e3,3)];  E(eg(6),:) = [v2,flag_E(e3,3)];
            ET(eg(5),1) = tri1; ET(eg(6),1) = tri3;
            flag_E(e3,1) = eg(5);  flag_E(e3,2) = eg(6); 
        else % this edge has been divided already
            eg(5) = flag_E(e3,2);  eg(6) = flag_E(e3,1);
            ET(eg(5),2) = tri1; ET(eg(6),2) = tri3;
        end
        eg(1) = posE + 1; eg(2) = posE + 2; posE = posE + 2;
        E(eg(1),:) = [flag_E(e2,3),flag_E(e3,3)];
        E(eg(2),:) = [flag_E(e3,3),v3];
        ET(eg(1),:) = [tri1;tri2]; ET(eg(2),:) = [tri2;tri3];
        ET(e1,find(ET(e1,:)==tri)) = tri3;
        % C: refine this triangle
        T(tri1,:) = [v1,flag_E(e3,3),flag_E(e2,3)];  TE(tri1,:) = [eg(1),eg(4),eg(5)];
        T(tri2,:) = [v3,flag_E(e2,3),flag_E(e3,3)];  TE(tri2,:) = [eg(1),eg(2),eg(3)];
        T(tri3,:) = [v2,v3,flag_E(e3,3)];  TE(tri3,:) = [eg(2),eg(6),e1];
    end
end

% run over the third category:
% scheme = 1;
for t = 1:length(T3)
    tri = T3(t);
    % find the longest edge as first edge:
    loc_1 = 1; longest = (V(T(tri,2),1)-V(T(tri,3),1))^2+(V(T(tri,2),2)-V(T(tri,3),2))^2;
    len = (V(T(tri,1),1)-V(T(tri,3),1))^2+(V(T(tri,1),2)-V(T(tri,3),2))^2;
    if len > longest
        loc_1 = 2; longest = len;
    end
    len = (V(T(tri,1),1)-V(T(tri,2),1))^2+(V(T(tri,1),2)-V(T(tri,2),2))^2;
    if len > longest
        loc_1 = 3; longest = len;
    end
    loc_2 = mod(loc_1,3)+1; loc_3 = mod(loc_2,3)+1;
    v1 = T(tri,loc_1); v2 = T(tri,loc_2); v3 = T(tri,loc_3);
    % A: add triangle first
    tri1 = posT + 1; tri2 = posT + 2;
    tri3 = posT + 3; tri4 = tri; 
    posT = posT+3; 
    % B: cope with edge
    e1 = TE(tri,loc_1); e2 = TE(tri,loc_2); 
    eg = zeros(9,1);   e3 = TE(tri,loc_3);
    if flag_E(e1,1) == 0 % this edge have not been divided yet, refine it!
        eg(1) = posE+1; eg(2) = e1; posE = posE + 1;
        E(eg(1),:) = [v2,flag_E(e1,3)];  E(eg(2),:) = [v3,flag_E(e1,3)];
        ET(eg(1),1) = tri1; ET(eg(2),1) = tri4;
        flag_E(e1,1) = eg(1);  flag_E(e1,2) = eg(2); 
    else % this edge has been divided already
        eg(1) = flag_E(e1,2);  eg(2) = flag_E(e1,1);
        ET(eg(1),2) = tri1; ET(eg(2),2) = tri4;
    end
    if flag_E(e2,1) == 0 % this edge have not been divided yet, refine it!
        eg(3) = posE+1; eg(4) = e2; posE = posE + 1;
        E(eg(3),:) = [v3,flag_E(e2,3)];  E(eg(4),:) = [v1,flag_E(e2,3)];
        ET(eg(3),1) = tri4; ET(eg(4),1) = tri3;
        flag_E(e2,1) = eg(3);  flag_E(e2,2) = eg(4); 
    else % this edge has been divided already
        eg(3) = flag_E(e2,2);  eg(4) = flag_E(e2,1);
        ET(eg(3),2) = tri4; ET(eg(4),2) = tri3;
    end
    if flag_E(e3,1) == 0 % this edge have not been divided yet, refine it!
        eg(5) = posE+1; eg(6) = e3; posE = posE + 1;
        E(eg(5),:) = [v1,flag_E(e3,3)];  E(eg(6),:) = [v2,flag_E(e3,3)];
        ET(eg(5),1) = tri2; ET(eg(6),1) = tri1;
        flag_E(e3,1) = eg(5);  flag_E(e3,2) = eg(6); 
    else % this edge has been divided already
        eg(5) = flag_E(e3,2);  eg(6) = flag_E(e3,1);
        ET(eg(5),2) = tri2; ET(eg(6),2) = tri1;
    end
    eg(7) = posE + 1; E(eg(7),:) = [flag_E(e1,3),flag_E(e2,3)];
    eg(8) = posE + 2; E(eg(8),:) = [flag_E(e1,3),v1];
    eg(9) = posE + 3; E(eg(9),:) = [flag_E(e1,3),flag_E(e3,3)];
    posE = posE + 3;         ET(eg(7),:) = [tri3,tri4];
    ET(eg(8),:) = [tri2,tri3]; ET(eg(9),:) = [tri1,tri2];
    % C: refine this triangle
    T(tri1,:) = [flag_E(e3,3),v2,flag_E(e1,3)]; T(tri2,:) = [flag_E(e3,3),flag_E(e1,3),v1];
    T(tri3,:) = [flag_E(e2,3),v1,flag_E(e1,3)]; T(tri4,:) = [flag_E(e2,3),flag_E(e1,3),v3];
    TE(tri1,:) = [eg(1),eg(9),eg(6)]; TE(tri2,:) = [eg(8),eg(5),eg(9)];
    TE(tri3,:) = [eg(8),eg(7),eg(4)]; TE(tri4,:) = [eg(2),eg(3),eg(7)];
end      
%             % At last update the information of areas:
%             x21 = V(T(1:posT,2),1) - V(T(1:posT,1),1);
%             x31 = V(T(1:posT,3),1) - V(T(1:posT,1),1);
%             y21 = V(T(1:posT,2),2) - V(T(1:posT,1),2);
%             y31 = V(T(1:posT,3),2) - V(T(1:posT,1),2);
%             areas = (x21.*y31 - x31.*y21)/2;